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In the context of subaqueous ripple and dune formation, we present here a Reynolds aver- 
aged calculation of the turbulent flow over a topography. Using a Fourier decomposition 
of the bottom elevation profile, we perform a weakly non- linear expansion of the veloc- 
ity field, sufficiently accurate to recover the separation of streamlines and the formation 
of a recirculation bubble above the some aspect ratio. The normal and tangential basal 
stresses arc investigated in details; in particular, we show that the phase shift of the shear 
stress with respect to the topography, responsible for the formation of bedforms, appears 
in an inner boundary layer where shear stress and pressure gradients balance. We study 
the sensitivity of the calculation with respect to (i) the choice of the turbulence closure, 
(ii) the motion of the bottom (growth or propagation), (iii) the physics at work in the 
surface layer, responsible for the hydrodynamic roughness of the bottom, (iv) the aspect 
ratio of the bedform and (v) the effect of the free surface, which can be interpreted in 
terms of standing gravity waves excited by topography. The most important effects are 
those of points (iii) to (v), in relation to the intermixing of the different length scales of 
the problem. We show that the dynamical mechanisms controlling the hydrodynamical 
roughness (mixing due to roughness elements, viscosity, sediment transport, etc) have an 
influence on the basal shear stress when the thickness of the surface layer is comparable 
to that of the inner layer. We evidence that non-linear effects tend to oppose linear ones 
and are of the same order for bedform aspect ratios of the order of 1/10. We show that 
the influence of the free surface on the basal shear stress is dominant in two ranges of 
wavelength: when the wavelength is large compared to the flow depth, so that the inner 
layer extends throughout the flow and in the resonant conditions, when the downstream 
material velocity balances the upstream wave propagation. 



1. Introduction 

The formation of ripples and dunes at the surface of an erodiblc sand bed results from 
the interplay between the relief, the flow and the sediment transport. The aim of these 
two companion papers is to propose a coherent and detailed picture of this phenomenon 
in the generic and important case of a unidirectional turbulent stream. This first part is 
devoted to the study of the stationary flow over a wavy rough bottom. In the second part 
we propose a common theoretical framework for the description of the different modes 
of sediment transport. Hydrodynamics and transport issues at hand, we then revisit the 
linear instability of a flat sand bed submitted to a water shear flow and show that, in 
contrast to ripples, subaqueous dunes cannot form by a primary linear instability. 
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Figure 1. Iso- velocity lines over a wavy bottom (data after |Poggi et al. 20 07). The fluid 
is flowing from left to right along the a;-axis. z is perpendicular. The bottom profile is 
2 = Z(x) = £cos(fcc). In that experiment, the wavelength and amplitude of the bumps are 
\ — 2n/k — 3.2 m and £ = 0.08 m, for a water depth H = 0.6 m). The point of maximum shear 
on the bump (where the lines are squeezed) is located upstream the crest. 



It has long been recognised that the mechanism responsible for the formation and 
growth of bedforms is related to the phase-lag between sediment transport and bed eleva- 
tion ( |Kennedy 1963||Reynolds 1965||Kennedy 1969| [ Smith 1970l[Hayashi 1970|IParker 19751 
Engclund & Fredsoc 1982 McLean 1990]) . It has been shown in the context of aeolian 
dunes that this lag comes from two contributions, which can be considered as independent 
as the time scale involved in the bed evolution is much slower than the hydrodynamics 
relaxation (lAndreotti et al. 2002|.|Kroy et al. 2002[ IValance 2005|) . First there is a shift 
between the bed and the basal shear stress profiles. This shift purely results from the 
hydrodynamics and its sign is not obvious a priori, i.e. the stress maximum can be 
either upstream or downstream the bed crest depending on the topography or the prox- 
imity of the free surface. The second contribution comes from the sediment transport: 
the sediment flux needs some time/length to adapt to some imposed shearing. This 
relaxation mechanism induces a downstream lag of the flux with respect to the shear. 
When the sum of these two contributions results in a maximum flux upstream the bed 
crest, sediment deposition occurs on the bump, leading to an unstable situation and thus 
to the amplification of the disturbance. In part 1, we shall focus on the first of these 
contributions, the second one being treated in part 2. 

We consider here the generic case of a flow over a fixed sinusoidal bottom of wave- 
length A (see figure [T] for an illustration of the geometry and some notations). In order 
to obtain the basal shear stress and in particular its phase shift with respect to the to- 
pography, the equations of hydrodynamics must be solved in this geometry. The case of 
viscous flows has been investigated by Benjam in 1959 Bordn er~19781|Caponi et al~~l 982, 
ICharru fc Hinch 20001 Lagree 2003HValance fc La nglois 2005} The first attempts to model 
the high Reynolds number regime in the context of ripples and dunes in rivers have dealt 
with potential flows ( |Kennedy 1963[ [Reynolds 1965[IColeman fc Fenton 2000[) . for which 
the velocity field does not present any lag with respect to the bottom. The shallow-water 
approximation (Gradowczyk 1970) implies that the bedforms spread their influence on 
the whole depth of the flow. However, patterns only have a significant influence within a 
vertical distance on the order of their wavelength. It is then crucial to compute explicitly 
the vertical flow structure, taking into account the turbulent fluctuations. 

In order to overcome the flaws of the perfect flow, constant eddy viscosity closures have 
been tried to improve Kennedy's original model (Engelund 1970, Smith 1970, Frc ds0e 1974| ). 
Further progress has been made by IRich ards 1980] who used a more sophisticated mod- 
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clling with an additional equation on the turbulent energy and a closure which involves 
a Prandtl mixing length in the expression of the eddy viscosity. |Sumer fcTB akioglu 1984 
made use of the same turbulent modelling, but in the case of an infinite water depth. 
A mixing length approach was also used by[Kobayashi & Madscn 1985 to improve Ben- 
jamin's laminar description. 



In the meteorological context of atmospheric flows over low hills, a deep and fun- 
damental understanding of the physics of turbulent flows over a relief has been devel- 
oped from the 70's (see the review bv IBelcher fc Hunt 1998|) . Starting with the semi- 
nal work of I Jackson fc Hunt 1975[ further refined by|Sykcs 1980 and IHunt et al. 1988|, 
the gross emerging picture is that the flow can be thought of as composed of two (or 
more) layers, associated with different physical mechanisms and different length scales. 
Jac kson fc Hunt 19751 have been able to compute analytically the basal shear stress for 
asymptotically large patterns, under an infinite flow depth assumption. Their ideas 
have been discussed in a rather vast literature. The predictions of these calculations, 
and in particular this layered structure of the flow, has been compared with exper- 
iments (see e.g. IBritter et al. 19811 |Gong fc Ibbetson 1989[ |Finnigan et al. 1990| , or 
field measurements on large scale hills (see e.g. the review paper by Taylor et al. 1987] ), 
with a good degree of success, especially on the upstream side of the bumps. More- 
over, they have been tested against the results of the numerical integration, in vari- 
ous configurations, of Navier-Stokes equations closed with different turbulent closures 
dTaylor 1977a] |Taylor 1977b] |Richards fc Taylor 1981j |Ayotte et al. 1994| ). The rele- 
vance of this approach for the description of the flow and the stresses around aeolian 
sand dunes has also been investigated (see e.g. |Weng et al. 199lj ), and is amongst the 
current directions of research in that community (Wiggs 2001). 



Because the prediction of the stable or unstable character of a flat sand bed submitted 
to a turbulent shear flow is very sensitive to the way both hydrodynamics and transport 
issues are described and intermixed, we find it useful to discuss at length, in the two 
parts of this paper, the different mechanisms and scales involved at the different steps of 
the modelling. It is indeed particularly revealing that, despite the fact that the approach 
presented here is very close to those of Richards 1980 and Colombini 2004 and has been 
motivated by these works, we basically disagree with their conclusions, especially that 
river dunes are initiated by the linear instability of a flat bed. The detail discussion 
we provide here gives also the opportunity to revisit the still debated question of the 
subaqueous ripple size selection (see ICharru 20061 and references therein) , as well as the 
important issue of the classification of bedforms ( [Ashley 1990p . 



This article is structured as follows. In the next section, we briefly recall the equations 
for the base flow over a uniform bottom. We then study the linear solution in the case 
of wavelengths much smaller than the flow depth. Importantly, the sensitivity of these 
linear results with respect to various changes in the modelling is tested in sections [4] and 
OH Section [6] is devoted to the derivation of the first non-linear corrections. In section 
[7] we investigate the effect of the free surface in the case of wavelengths comparable 
or larger than the flow depth and interpret it in terms of topography induced standing 
gravity waves. Finally, we provide in section [8] a qualitative summary of the main results 
of the paper. The most technical considerations are gathered in appendices. 
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2. Turbulent flow over a uniform bottom 

2.1. The logarithmic law 

We consider a turbulent flow over a relief. Following Reynolds' decomposition between 
average and fluctuating (denoted with a prime) quantities, the equations governing the 
mean velocity field itj can be written as: 

d iUi = 0, (2.1) 
D t Ui = d t Ui + u-jdjUi = -djT i:i - dtp, (2.2) 



where = u^u'- is the Reynolds stress tensor ( |Reynolds 1874| ). For the sake of simplicity, 
we omit the density factor p in front of the pressure p and the stress tensor. The aim of 
this paper is to describe quantitatively the average flow over a fixed corrugated boundary 
within this framework. The reference state is the homogeneous and steady flow over a flat 
bottom, submitted to an imposed constant shear stress t xz = —u^. The turbulent regime 
is characterised by the absence of any intrinsic length and time scales. At a sufficiently 
large distance z from the ground, the only length-scale limiting the size of turbulent 
eddies - the so-called mixing length L - is precisely z; the only mixing time-scale is given 
by the velocity gradient |9 z tta;|. As originally shown by Prandtl 1925 , it results from this 
dimensional analysis that the only way to construct a diffusive flux is a turbulent closure 
of the form: 

Txz = -K 2 L 2 \d z u x \d z u x , (2.3) 

where the mixing length is L — z and k ~ 0.4 is the (phenomenological) von Karman 
constant. After integration, one obtains that the velocity has a single non zero component 
along the x-axis, which increases logarithmically with z (Tritton 1988): 

^ = -ln(-V (2.4) 

where zq is a constant of integration called the hydrodynamical roughness. This expres- 
sion does not apply for z — > 0. There should be layer of thickness h close to the bottom, 
called the surface layer, matching the logarithmic profile to a null velocity on the ground. 

2.2. Hydrodynamical roughness 

The hydrodynamical roughness zq should be distinguished from the geometrical (or phys- 
ical) roughness of the ground, usually defined as the root mean square of the height 
profile variations. z is defined as the height at which the velocity would vanish, when 
extrapolating the logarithmic profile to small z. The physical mechanism controlling 
zq can be of different natures. If the ground is smooth enough, a viscous sub-layer of 
typical size 0(v/u*) must exist, whose matching with the logarithmic profile determines 
the value of zq. On the contrary, if the geometrical roughness is larger than the vis- 
cous sub-layer, turbulent mixing dominates at small z with a mixing length controlled 
by the ground topography. In the case of a static granular bed composed of grains 
of size d, reported values of the hydrodynamical roughness are reasonably consistent 
(zo ~ rf/30 in |Bagnold 1941[ z ~ d/24 in |Schlichting fc Gersten 2000| and z ~ d/10 



dKamphuis 1974 I Andreotti 2004jl . In section |6| we will justify the connection between 



geometrical and hydrodynamical roughness on a rigourous basis and show that they are 
not simply proportional. 

The situation is of course different in the presence of sediment transport, which may 
(or not) induce some negative feedback on the flow. In this case, the hydrodynamical 
roughness z Q may directly be controlled by the transport characteristics (e.g. mass flux 
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and grain trajectories). Nature presents many other physical processes controlling the 
roughness: for instance, the flexible stems of wetland plants in low marshes or, for the 
wind, the canopy or the waves over the ocean. In all these cases, it can be assumed 
that the logarithmic law is a good approximation of the velocity profile above the surface 
layer, with a single known parameter zq. 

We will first consider the asymptotic limit in which the typical relief length - say, the 
dune wavelength A - is much larger than the surface layer thickness h$. The relief is 
locally flat at the scale ho, so that there must be a region close to the ground where 
the velocity profile shows a logarithmic vertical profile. We will then discuss the case of 
moderate values of the ratio X/Hq, for which the flow becomes sensitive to the details of 
the mechanisms controlling the roughness. 

2.3. A turbulent closure 
In the logarithmic boundary layer, the normal stresses can be written as: 

_ 1 

z ~ 3 



„,, - , . ; - -Tu with t u = K 2 x 2 L 2 \d z u x \ 2 , (2.5) 



where x is a second phenomenological constant estimated in the range 2.5 — 3. Note that 
X does not have any influence on the results as it describes the isotropic component of 
the Reynolds stress tensor, which can be absorbed into the pressure terms. Normal stress 
anisotropy is considered in section [5] and appendix [X] Introducing the strain rate tensor 
jij = djUj + dj Uj and its squared modulus I7I 2 = -^jijjij, we can write both expressions 



(2.3 1 and (2.5 1 in a general tensorial form: 



T ii = K 2 L*\M\±x'\i\5 ij -i ij \. (2.6) 

In this paper, we focus on 2D steady situations, i.e. on geometries invariant along 
the transverse y-direction, see figure [l] As they are of permanent use for the rest of 
the paper, we express the components of the velocity and stress equations in the x- and 
z-directions. The Navier-Stokes equations read: 

d x u x + d z u z = 0, (2.7) 
u x d x u x + u z d z u x = -d x p - 8 z t xz - d x T xx , (2.8) 
u x d x u z + u z d z u z = -d z p - 3 z t zz - d x T zx . (2.9) 

The stress expressions are the following: 

t xz = -n 2 L 2 \j\i xz , (2.10) 

t xx = -K 2 L 2 m xx + i K 2 x 2 L 2 |7| 2 , (2.11) 

T zz = -n 2 L 2 m zz + l r 2 X 2 L 2 \^\ 2 . (2.12) 

In these expressions, the strain tensor components are given by 

jxz = izx = d z u x + d x u z , j xx = 2d x u x and A f zz = 2d z u z = -j xx , (2.13) 
and the strain modulus by: 

| 7 | 2 = 2{d x u x f + 2(0 z u z ) 2 + (d z u x + 8 x u z f = A{8 x u x ) 2 + (d z u x + 8 x u z ) 2 . (2.14) 
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3. Unbounded turbulent boundary layer over a wavy bottom 

We now consider the turbulent flow over a wavy bottom constituting the floor of an 
unbounded boundary layer. In rivers, this corresponds to the limit of a flow depth H 
much larger than the bed- form wavelength A. The solution is computed as a first order 
linear correction to the flow over a uniform bottom, using the first order turbulent closure 
previously introduced. 



3.1. Linearised equations 
For small enough amplitudes, we can consider a bottom profile of the form 

Z(x) = (cos(kx) 



(3.1) 



without loss of generality. A = 27r/fc is the wavelength of the bottom and C the amplitude 
of the corrugation, see figure [T| The case of an arbitrary relief can be deduced by a 
simple superposition of Fourier modes. We introduce the dimcnsionless variable r\ = kz, 
the dimensionless roughness ?/ = kz and the function: 



In 



(3.2) 



We also switch to the standard complex number notation: Z(x) = (^e %kx (real parts of 
expressions are understood). 



We wish to perform the linear expansion of equations ( 2.7 )-( 2.14) with respect to the 



small parameter k(. The mixing length is still defined as the geometrical distance to the 
bottom: L = z — Z. We introduce the following notations for the two first orders: 

u, [n + k(e lkx U] , (3.3) 

u*k(e ikx W, (3.4) 

t zx = -u\ [1 + k(e ikx S t ] , (3.5) 
'1 



P 



PO 



kQe lKX S t \ 
- k(e lkx S„ 



kC,e lkx S. 



ikx Q 



(3.6) 
(3.7) 
(3.8) 



The quantities U, W, etc, are implicitly considered as functions of r\. An alternative 
choice is to consider functions of the coordinate £ = rj — kZ. Such alternative functions are 
denoted with a tilde to make the distinction. This important - but somehow technical - 
issue of the choice of a representation is discussed in appendix[C] Although the curvilinear 
and Cartesian systems of coordinates are equivalent, the distinction between the two is 
of importance when it comes to the expression of the boundary conditions, and for the 
range of amplitudes £ for which the linear analysis is no more valid (see section |6| . In 
particular, vertical profiles in the forthcoming figures will be mostly plotted as a function 
of the shifted variable £. 

The linearised strain rate tensor reads 



%z - j zx = ku«n' + u«k 2 (e lkx (U' + iW), 
j xx = 2iu.,k 2 (e lkx U 1 
% z = 2u*k 2 (e lkx W, 

\M = \ixzl 



(3.9) 
(3.10) 
(3.11) 
(3.12) 
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and the stress equations can be simplified into 

fjfSt = 2(U' + iW) - 2/t 2 r?// 3 , 
(j,'S xx = ~2iU + 2 -x\U' + iW) - 2 - X 2 ^'\ 

»'S ZZ = -2W' + 2 - X 2 {U' + iW) - 2 - X 2 ^' 2 . 
Finally the Navier-Stokes equations lead to 

W' = -iU, 

S' t = fiiU + fi'W + iS n + iS xx - iS zz , 
S' n = -niW + iSt. 



Taking the difference of equations (3.141 and (3.151, one can compute 



Sxx S zz — 



-4iU 



to obtain four closed equations: 

U' = -iW+^fjfSt + kh' 2 , 
W = -iU, 

S' t = +j}jU + li'W + iS n , 
S' n = -ipW + iSt. 



(3.13) 
(3.14) 

(3.15) 



(3.16) 
(3.17) 
(3.18) 



(3.19) 

(3.20) 
(3.21) 
(3.22) 
(3.23) 



Introducing the vector X = (U,W, St, S n ), we finally get at the first order in k( the 
following compact form of the equation to integrate: 



/ -i \p! 0\ 



di] 



X = VX + S, with V 



and S 



(3.24) 





W+p) i 

\ -ni i OJ 

The general solution of this equation is the linear superposition of all solutions of the 
homogeneous system (i.e. with 5 = 0), and a particular solution X a . 






V o J 



3.2. Boundary conditions 



Four boundary conditions must be specified to solve the above equation (3.241 
upper boundary corresponds to the limit rj 



The 



oo, for which we ask that the vertical 
fluxes of matter and momentum vanish asymptotically. This means that the first order 
corrections to the shear stress and to the vertical velocity must tend to zero: W(co) = 
and 5 t (oo) = 0. In practice, a boundary at finite height H (at j]h = kH) is introduced, 
at which we impose a null vertical velocity W(t]h) = and a constant tangential stress 
pu 2 so that St(r)H) — 0. This corresponds to a physical situation where the fluid is 
entrained by a moving upper plate, for instance a stress-controlled Couette annular cell. 
Then, we consider the limit H — > +oo, i.e. when the results become independent of H . 

The lower boundary condition must be specified on the floor (77 — > kZ). We consider 
here the limit in which the surface layer thickness ho is much smaller than the wavelength 
A. This allows to perform an asymptotic matching between the solution and the surface 
layer, whatever the dynamical mechanisms responsible for the hydrodynamical roughness 
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Zq are. Indeed, focusing on the surface layer, we know that in the limit z S> ho, the 
asymptotic behaviour of the local tangential velocity u should be a logarithmic profile 



controlled by the local shear stress r and the roughness zO. The solution of (3.241 should 
thus match this asymptotic behaviour as 77 — > kZ. Thus, zq is the only parameter 
inherited from the surface layer in the limit ho <C A. We will investigate the situation 
where this approximation is not valid anymore in sec tion [4] 

In the limit z—Z <C A, the homogeneous solution of (3.24 ) can be expanded in powers of 



77 and In and expressed as the sum over four modes. Adding the asymptotic behaviour 



of the particular solution X s = (- 



X 



KTj ' 



■no 



0,0), the full solution writes: 



/ M 2 /4 \ 




( M/2 ^ 




( 1 ^ 




^ 2 /(4«) 


1 


+ a 2 


—i 77 /Lt/2 




— i 77 




A* 


1 


\-v 2 Kv) J 


i 7] 






\ in / 






{ 1 / 



-)• 



(3.25) 



The next terms in this expansion are ©(fyln 

The values of the four coefficients a\, 04 are selected by the matching with the 
surface layer. a\ would correspond to a non vanishing normal velocity through the 
surface layer and should thus be null. a 2 precisely corresponds to the logarithmic profile 
with a roughness Zq and a basal shear stress modulation a 2 . This gives a 2 = S't(O). 0,3 
would correspond to a modulation of the local roughness - more precisely of its logarithm. 
We do not consider such a modulation so that 03 = 0. 04 corresponds to a sub-dominant 
behaviour associated to the basal pressure modulation (04 — S n (Q)). In summary, the 
functions U, W, St and S n should follow the following asymptotic behaviour: 



U{rj) 



St(0) 
2k 



In 



iS n (0) 



2 k 



->1 



w , , iS t (0) 
W(r)) = 



2k 
S n (0). 



iS, 



ln^ 

m 



- 1 



1 

KT] 

■Sn(O ), 

4k 



-V 



In 



no 



(3.26) 

(3.27) 

(3.28) 
(3.29) 



The region of thickness I in which this asymptotic behaviour constitutes a good ap- 
proximation of the flow field is called the inner layer. Equation (3.291 means that the 
total pressure p = p + t/j/3 is constant across this boundary layer: 



d z p = 



(3.30) 



and equation (3.281 that the shear stress decreases linearly with height according to 

d x p + d z T xz = 0, 



(3.31) 

The tangential pressure gradient is balanced by the normal shear stress, which means that 
inertial terms are negligible or equivalently that the fluid is in local equilibrium. In terms 
of energy, the space variation of the internal energy (pressure) is dissipated in turbulent 
" friction" . These two equations correspond to the standard lubrication approximation 
for quasi-parallel flows. 

3.3. Equations solving 

In practice, we solve the equations using a fourth order Runge-Kutta scheme with a 
logarithmic step. The integration is started at an initial value of 77 inside the inner layer 
i.e. which verifies 77 In 2 ^ < 1). We write the solution as a linear superposition of the 
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Figure 2. Vertical profiles of the first order corrections to velocities and stresses for r/o — 1CP 4 . 
f = rj — kZ is the distance to the bottom, rescaled by the wavenumber. In all panels, the solid 
lines represent the real parts of the functions, whereas the dotte d lines represent the imaginary 
ones. Dashed lines show the asymptotic behaviours ( 3.26)3.29 1 used as boundary conditions. 
They match the solutions in the inner layer, which extends up to r/ ~ k£ ~ 10~ 2 here. We note 
St(0) = A + iB and S n {0) = C + iD. Close to the boundary, a plateau of constant shear stress 
can be observed, which corresponds to the logarithmic zone. It is embedded into a slightly larger 
zone of constant pressure in which the shear stress varies linearly. 



form X = X s + S t (0)X t + S n (0)X n , where the different terms verify: 



ch] 



X s = VX S + S starting from ^(77) 



( \ 

K r/o 


V ) 



(3.32) 
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Figure 3. Flow streamlines above a wavy bottom of rescaled amplitude fc£ = 0.5 (aspect ratio 
~ 1/6), computed from the linearised equations (770 = 10 -4 ). The flow direction is from left to 
right. Note the left-right asymmetry of the streamlines around the bump in the inner layer (grey 
lines). Note also the onset of emergence of a recirculation bubble in the troughs. The thick line 
in the top right corner shows the positions that maximises the velocity along a streamline. 



drj 



X t = VX t starting from X t (77) = 



2k 770 

a (In -a 

« V Vo 



- 1 



drj 



X n = VX n starting from X n (rj) 



3L 
4 k. 

V 

V 1 J 



(3.33) 



(3.34) 



The boundary conditions on the bottom arc then automatically satisfied, and the top 
ones give algebraic equations on the real and imaginary parts of St(0) and S n (0), which 
can be solved easily. We have checked that the result is independent of the initial value 
of f], as long as it remains in the announced range. 



3.4. Results 



The velocity and stress profiles resulting from the integration of equation (3.24) are 
displayed in figure [2] Looking at panel (c), one can clearly see the region close to 
the bottom where the shear stress is constant, while the horizontal velocity component 
(panel a) exhibits a logarithmic behaviour. This plateau almost coincides with the inner 
layer, which is the zone where the solution is well approximated by the asymptotic 
behaviour derived above. The inner layer is embedded in a wider region characterised by 
a constant pressure (panel d). The estimate of the thickness i is of crucial importance 
for the transport issue (see next section and Part 2). t is the scale at which inertial terms 
are of the same order as stress ones in the Reynolds averaged Navier-Stokes equations. 
The original estimation of i given by [Jackson & Hu nt 19751 was further discussed in 
several later papers (see e.g. |Taylor et al. 1987] [Clausscn 1988, Belja ars~fe Taylor~ 989, 
Finnigan et al. 19"90|. Our data are in good agreement with the scaling proposed by 



Taylor et al. 1987 



1 1 \ n 2 - = 0{l). 



(3.35) 



A zo 

Consistently, this scaling relationship is precisely that of the first neglected terms in the 
asymptotic expansion (3.251. Away from the bottom, all profiles tend to zero, so that 
one recovers the undisturbed flow field (2.4 1 at large 77. The shape of these profiles are 



very consistent with the work of |Ayotte et al. 1994[ who have compared the influence of 
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u u 

Figure 4. Vertical profiles of the first order correction to the horizontal velocity for 770 = 10~ 4 . 
(a) Lin-lin plot in the shifted representation, (b) Lin-Log plot in the non-shifted representation. 
The solid lines correspond to the real part and the dotted line to the imaginary one. The velocity 
disturbance decreases exponentially over one wavelength (dashed line) . In this outer region, the 
Reynolds stress can be neglected. 



the closure scheme on the linear flow over a relief, which means that the precise choice 
of the turbulent closure does not affect significantly the results. 

In order to visualise the effect of the bottom corrugation on the flow, the flow stream- 
lines are displayed in figure[3](see appendix [D] for explanations about their computation). 
It can be observed that the velocity gradient is larger on the crest than in the troughs 
as the streamlines are closer to each other. The flow is disturbed over a vertical distance 
comparable to the wavelength. A subtler piece of information concerns the position along 
each streamline at which the velocity is maximum. These points are displayed in the right 
corner of figure [3] Away from the bottom, they are aligned above the crest of the bump. 
Very close to it, however, they are shifted upstream. In other words, the fluid velocity is 
in phase with the topography in the upper part of the flow, but is phase advanced in the 
inner boundary layer where the shear stress tends to its basal value. In this inner layer, 



the profile is well approximated by its asymptotic expression (3.26). 

An inspection of the velocity profile evidences two distinct regions (see figures [2] and [4]), 
in which we recognise those at the basic partitioning of the flow in Jackson & Hunt work 
(1975), and subsequent papers. There is an outer region (77 ^> k€), where U decreases 
exponentially with rj (figure ptb)). Seeking for asymptotic solutions decreasing as e~ CT? ', 
one has to solve the eigenvalue problem VX = —oX for asymptotically large values of 
rj. At the two leading orders, the decrease rate a is given by: 

2i (a 4 + 1) K 2 r/ + (a 2 - l) In — = 0. (3.36) 

m 

The asymptotic behaviour is an oscillatory relaxation corresponding to a = (1 ± i)/\/2. 
However, the observed decrease corresponds to the intermediate asymptotic regime rj < 
In — for which the solution is a = 1. This behaviour is reminiscent from that of an 

Vo 

inviscid potential flow. In other words, the effect of the turbulent shear stress on the 
flow disturbance can be neglected. 

The intermediate region between the inner and the outer layers is responsible for the 
asymmetry of the flow as well as the upstream shift of the maximum velocity discussed 
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Figure 5. Coefficients A, B, C and D as a function of rjo = kzQ. These piots show the depen- 
dence of the basal shear and normal stresses with the number of decades separating the wave- 
length A from the soil roughness zo, for a given bump aspect ratio. The solid line corresponds 
to the results of the model, using the asymptotic matching with the surface layer. The dashed 
lines represent the analytical formula deduced from I Jackson fc Hunt 19751 by [Kroy et al. 2002| 
They agree well at very small r/o. 



above. Let us emphasise again that this is the physical key point for the formation of 
bcdforms. One can understand the reason of the phase shift with the following argument. 
The external layer can be described as a perfect irrotational flow. Since the elevation 
profile is symmetric, the streamlines are symmetric too, as the flow is solely controlled 
by the balance between inertia and the pressure gradient induced by the presence of the 
bump. As a consequence, the velocity is maximum at the vertical of the crest. Now, 
inside the inner layer, this flow is slowed down by turbulent diffusion of momentum. 
Focusing on the region of matching between these outer and inner regions, the velocity 
needs some time to re-adapt to a change of shear stress, due to inertia. Thus, the shear 
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stress is always phase-advanced with respect to the velocity. One concludes that the 
basal shear stress is phase-advanced with respect to the bump. 

As mentioned in the introduction, we are especially interested in the shear stress and 
pressure distributions on the bottom. We note <St(0) = A + iB and S n (0) = C + iD. The 
ratio B/A is the tangent of the phase shift between the shear stress and the topography. 
It is positive as the shear stress is phase advanced. The four coefficients A, B, C and D 
are displayed as a function of rjo in figure [5] Their overall dependence with rjo is weak, 
meaning that the turbulent flow around an obstacle is mostly scale invariant. More 
precisely, following Jackson & Hunt's work (Jacks on fc Hunt 19751 |Kroy et al. 2002| ), it 
has been shown that, for asymptotically small 770, one expects logarithmic dependencies: 

In 2 ($ 2 /ln$) / vr \ ln 2 ($ 2 /ln$) 

A= \ ' > (l + ln^ + 21n- + 47 g ) and B = it \ ' 1 , (3.37) 

2 In <p v Z / 2 In <p 

where Euler's constant is je — 0.577, (j> is defined by the equation <j> In (j> = 2k 2 <I> and with 
$ = 7r/(2r/o). Note that A tends to 2 and B to as 770 — » 0, as expected when the inner 
layer thickness I vanishes. In this limit, the basal shear stress is directly proportional to 
the square of the velocity inherited from the outer layer, which is solution of the potential 
flow problem. 

These expressions agree well with our numerical results for very small 770. However, for 
realistic values of r/ , e.g. 1CP 4 < 77 < 1CP 2 , this approximation cannot be accurately 
used as it leads to errors of order one - note that Jackson & Hunt's expressions tend 
to diverge at larger 770 • In comparison to A and B, we observe that the normal stress 
coefficients C and D are more robust with respect to the details of the model. In the 
limit of a perfect flow, the pressure varies as the square of the velocity. Here, one 
needs to consider the velocity at the scale A of the perturbation, say u*fj>, where the 
logarithmic factor [i should be evaluated for 77 of order unity. From this argument, we 
predict that the pressure coefficient C should scale as the square of In 770 (a parabola in 
figure [5|, which is very accurately verified. More precisely, C — [/i(l/4)] 2 is an almost 
perfect approximation. Finally, it can be observed that the normal stress is also in phase 
advance with respect to the bottom profile. The coefficient D is positive and shows a 
linear variation with In t]q. 



4. Effect of the mechanisms controlling the hydrodynamical 
roughness 

So far, the computation of the velocity and stress fields has been obtained without 



any specification of the physics at the scale of Zq, as the integration of equation (3.241 
was started in the inner layer rather than on the bottom. This is of course possible only 
if this layer is sufficiently thick, i.e. if I (or A) is much larger than the thickness of the 
surface layer ho introduced in section [2j We now discuss several ways to describe the flow 
inside the surface layer, and investigate the subsequent effect on the shape of the stress 
coefficients as functions of 770- These coefficients should be independent of the physics at 
work in this surface layer when 770 is small enough, but we expect larger differences for 
larger values of 7? . 

We first present a convenient phenomenological model of geometrically induced rough- 
ness, which does not involve additional parameters. Because of its simplicity, it will be 
used in the next sections, as well as in the second part of the paper. The results will be 
compared to a rigourous treatment resulting from the weakly non-linear hydrodynamical 
calculation in section [6] We then consider the case of a viscous surface layer. Inspired 
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Figure 6. Shear stress coefficients A and B (das hed lines) computed with the phenomenological 
model of geometrically induced roughness (TOT). For comparison, the solid lines display the 
reference case shown in figure [5] 

from the aeolian transport properties, we finally discuss the focus point assumption as 
a possible way to describe the situation in which the surface layer is governed by the 
presence of sediment transport. 



4.1. Geometrically induced roughness 

For an hydrodynamically rough bottom, the "small scale" roughness elements are larger 
than the viscous sub-layer. They are submitted to a turbulent drag from the fluid and 
reciprocally, their presence slows down the flow. The exchanges of momentum in the 
surface layer are thus dominated by the turbulent fluctuations. Following Richards 1980 
and others, a convenient phenomenological model is to define the mixing length involved 

— Z . In this way, L is still essentially the 



in the turbulent closure (2.6 1 as L 



geometrical distance to the bottom, except that it cannot be smaller than the roughness 
length. This choice reflects in a intuitive manner the physical picture one can infer from 
experiments or simulations where square-shaped roughness elements are glued on a flat 
wall (see e.g. |Perry et al. 19 69). We will show in section ^ that this picture must 
be refined when dealing with blunt roughness elements of moderate aspect ratio - for 
instance the surface of a sand bed. 

With this expression for the mixing length, the integration of starting equations in the 
uniform and steady case gives 



u x = — In ( 1 



z 



(4.1) 



where the lower boundary condition u x = can now be taken in z — 0. This expression 



is well approximated by the pure logarithmic profile (eq. 2.4) as soon as z is larger than, 
say, few zq. In other words, for this model, Iiq ~ z . 

The above description of the linear analysis, and in particular the expression of the 



matrix V and the vector S involved in (3.241, in the case of a wavy bottom is still valid, 



but now with the following expression for the function 



Kv) = 1 In ( 1 + ^ 



(4.2) 



The solution of (3.24 1 can again be written as a linear superposition X = X s + S t (0)X t + 
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Figure 7. Shear stress coefficients A and B computed with a viscous surface layer. Dotted 
line: TZt = 125; Dashed line: 1Z t = 1. As in figure[6] the solid lines display the reference case of 
figure [5] 



S n (0)X n , where these three vectors are solutions of 



dr] 



X s = VX S + S with 



dr) 



X t = VX t with 



dr) 



X n = VX n with 



X.(0) = 



X t (0) = 



X n (0) = 



( K7] ^ 




V J 

(°\ 



1 




v 1 y 



(4.3) 



(4.4) 



(4.5) 



This decomposition ensures the requirement that both components of the velocity vanish 
on the bottom, leading to W(0) = and U(0) — —^'(0) = — 1/(kt]o). As in the previous 
section, the coefficients S t (0) and 5„(0) are found by the upper boundary conditions. 

The coefficients A and B resulting from this integration are displayed in figure [6] One 
can see that, for r/o < 10~ 3 , they are not very much different from those obtained in the 
previous section. However, one can notice significant differences for r]o > 10 -2 . As the 
mixing length in the surface layer is larger in this case (L ~ zq) than in the asymptotic 
case (L ~ z — Z), the turbulent 'diffusion' is more efficient. This results into a larger 
phase advance for the shear stress (Fig. |6]c). For practical purposes and for later use 
in the second part of this paper, a very good empirical fit of the coefficients A and B is 
obtained with 



A = 2 



ai + a 2 R + a 3 R 2 + G^i? 3 
1 + a 5 R 2 + aaR* 



and B = 



h + b 2 R + b 3 R 2 + b 4 R 3 
1 + b 5 R 2 + 6 6 i? 4 



(4.6) 



with {01,02,03,04,05,06} = {1.0702,0.093069,0.10838,0.024835,0.041603,0.0010625}, 
{&i, 6 2 , 63, 64, b 5 , b 6 } = {0.036989, 0.15765, 0.11518, 0.0020249, 0.0028725, 0.00053483} and 
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><— = u%, (4.7) 



4.2. A viscous surface layer 

In hydraulically smooth situations, it is natural to expect that, very close to the bottom, 
the flow must be laminar and thus described by the equation 

du 
dz 

whose solution is 

u 2 

u x (z) = ^z. (4.8) 
v 

We thus neglect here the possibility of a phase shift across the viscous surface layer. 
The transition from viscous to turbulent regime is governed by the Reynolds number 
and occurs at a typical value TZt — 125. The surface layer thickness can then be easily 
computed as = — y/7Z t . At z = ho, both viscous and turbulent expressions for the 
velocity must coincide: 

u h = Jfct = — In — . (4.9) 

K z 

From this equality, we can deduce the hydrodynamical roughness seen from the inner 
layer, due to this viscous surface layer: 

v 



z Q = — y/Kte-"^. (4.10) 

Us, 

In the case of a sand bed, the transition between the hydrodynamically smooth and 



rough regimes occurs when the viscosity induced roughness (eq. 4.101 is of the order of 
the geometrically induced roughness (zq ~ d/10). 



With the corresponding value for rj — kz , we solve equation (3.241 in the usual 
manner, writing the solution in the form of the linear superposition as described above, 
except that the integration is started at the initial value 7/ = kho, in which we impose 
that the velocity is parallel to the bed and equal to Uh- At linear order, this leads to 

U(kh Q ) = -n'(kh ), (4.11) 
W(kh ) = iu h /u* = ifj,(kh ). (4-12) 

The resulting shear stress coefficients A and B are displayed in figure [7J As one can 
expect, in comparison to the reference case, they are smaller for larger values of TZt, 
and all different curves collapse as 770 — > 0. The viscous diffusion of momentum is less 
efficient than that induced by turbulent fluctuations. Moreover, in the Stokes regime, 
for Reynolds numbers much smaller than 1, the kinematic reversibility leads to a shear 
stress in phase with the topography. Consistently, it can be observed in figure [7jc) that 
the phase advance is reduced in the hydrodynamically smooth regime. We will show in 
the second part of this article that this explains the fact that the wavelength at which 
ripples appear is larger as the Reynolds number decreases. 

Experiments in the hydraulically smooth regime have been performed bv lZilker et al. 19771 
Zilker fc Hanratty 1977[ Abrams fc Hanratty 1985} who measured the ionic mobility be- 



tween two nearby electrodes. This current is assumed to be related, without any spatial 
or temporal lag, to the basal shear stress. The measured phase shift between the signal 
and the bottom topography could reach values as high as 80°. this would correspond to 
B/A = tan(807r/180) ~ 5.67. Within the present model, the phase shift remains much 
lower than the measured 80°. This unexpected value has been interpreted as the signa- 
ture of a lag of the laminar-turbulent transition with respect to the Reynolds number 
criterion TZ — IZf Further experiments based on a different measure principle are needed 
to understand this discrepancy. 
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Figure 8. (a) Shear stress coefficients A and B computed in the presence of a 'focus point' at 
height /iq, where the velocity is u^, as a function of r/o • A and B are larger for larger values 
of Uh/u, (1, 2, 5 and 10). However, the ratio B/A is less sensitive to this parameter, up to 
7/0 — 10 -3 . Again, the solid lines display the reference case of figure [H] 



This viscous surface layer model is an effective way to take bedload transport into 
account. As a matter of fact, anticipating on the part 2 of this paper where the dynamical 
mechanisms governing the sediment transport are discussed, transported particles are not 
passive and exert a stress on the fluid. Close to the transport threshold, their influence on 
the flow is negligible. However, as their density increases, transport induces a negative 
feedback on the flow, which should be taken into consideration in the hydrodynamics 
description. The simplest model of multi-layer sheet flow would be a Newtonian fluid 
whose viscosity increases with the concentration of moving sediments. In this large shear 
velocity regime, one thus expects a decrease of the phase-lag responsible for the ripples 
instability and possibly, a restabilisation of the bed. 

4.3. The focus point assumption 

An alternative manner to take the feedback of the transport on the flow into account 
can be achieved in analogy with the aeolian case, which provides the archetype of such a 
situation. In this case, it has been shown that the moving grains slow down the flow in 
the transport layer, whose thickness ho is independent of the shear velocity u, . Note that 
in the subaqueous case, the transport layer thickness is observed to gently increase with 
the shear stress (jAbbott fc Francis 19 77 1 [Fernandez Luque fc v an Beck 1976| ) close to the 
threshold, in the erosion limited regime (see part 2). Above ho, the effect of the particles 
on the flow is negligible and one recovers the undisturbed logarithmic velocity profile, 
but with a roughness larger than that without transport. Below h 0l the flow velocity 
is reduced and is independent of u* flUngar fc Haff 19871 lAndreotti 2004j) . As shown 
experimentally by |Bagnold 194l[ the velocity vertical profiles measured for different shear 
velocities thus cross at the 'focus point' z = ho and u x = Uh- At this point we have 

^ = lln^, (4.13) 

which means that the effective roughness in the logarithmic region, due to this transport 
layer, is 

zo = h e~ KUh/u - (4.14) 
To determine the flow field in such a situation, the crucial point is to compare ho 



18 



A. Fourriere, P. Claudin and B. Andreotti 



with the thickness of the inner layer £, i.e. the size of the constant stress plateau (see 
figure [2]) . If ho is larger than £, it means that one cannot reduce the transport issue to a 
relationship between the sediment flux and the basal shear stress only. In that case, the 
whole vertical velocity profile, which depends on the entire bottom elevation, is involved. 
Conversely, for ho < £, one can account for transport by modifying the bottom boundary 
conditions as follows. Following what we have done in the previous sub-section, we can 
impose that the fluid velocity at z = Z + h is parallel to the bed and equal to . At 
the linear order, we then get: 

U(kh ) = -//(fc/io) (4.15) 
W(kh ) = iu h /u* = in(kh ). (4-16) 

The result of this choice is shown for the stress coefficients in figure [8] for various values 
of Uh/us,. A and B are larger for increasing focus velocities, or equivalently larger focus 
altitude. As in the viscous surface layer case, all curves collapse for 770 — > because I gets 



larger in this limit (see equation (3.35 1). Interestingly, as far as bedforms are concerned 
(see Part 2), the ratio B/A is much less sensitive to variations of Uh/u*, at least in the 
region 770 < 1CP 3 . 

Finally, it should be noted that the focus point model only applies to the momen- 
tum limited transport regime (see part 2). Close to the transport threshold, in the 
erosion limited regime, the feedback of the particle transport on the flow is negligi- 
ble and the transport should not to be taken into account in the hydrodynamical cal- 
culation. This approach should be distinguished from that used by IColombim~2 004 , 
Colombini & Stocchi no 20051 [Colombini & Stocchi no 20081 In those articles, the flow 
boundary conditions are taken on the bottom, below the transport layer (vanishing ve- 
locities), meaning that the feedback of the transport on hydrodynamics is neglected. 
However, the stresses are evaluated in ho, above the transport layer. This does not con- 
stitute a self-consistent choice. Moreover, this introduces a free parameter in the model 
which can be tuned to choose the phase shift at will. We will present in the second part 
experimental evidences that this choice is not correct. 

4.4. Concluding remarks 

For these three dynamical mechanisms controlling the hydrodynamical roughness zq, we 
have seen that the asymptotic regime is recovered when the surface layer thickness h 
is smaller than the inner layer thickness £. As I is much smaller than the wavelength 
A (for standard bedforms, lj\ = O(10 -3 )), this constitutes a rather restrictive condi- 
tion. Whenever ho is larger than £, specific hydrodynamic models should be derived to 
determine the relations between stresses and topography. 

The phase shift between the basal shear stress and the topography originates from 
the interplay between inertia and shear stress. The different models of surface layer 
correspond to different ways of mixing momentum in the direction normal to the wall. 
Although the argument is general, one sees that the precise value of the phase shift is 
rather sensitive to the physical origin of the momentum fluxes. In particular, viscous 
diffusion leads to a much smaller phase advance than turbulent mixing. 



5. Robustness of the results 

In the same spirit as the previous section, where we have investigated the effect of 
different ways to treat the surface layer on the stress coefficients A and B, we would 
like now to show the robustness of the results when considering a second order turbulent 
closure, Reynolds stress anisotropy or a moving bottom. 
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Figure 9. Effect of the parameter /3,which is a non dimensional parameter encoding the time 
lag between a change in the strain rate and that of the Reynolds stress, (a-c) Vertical profiles 
of the horizontal component of the velocity U = U + fi' for 770 = 10~ 4 for (a) (3 — 0, (b) (3 = 1 
and (c) (3 = 10 respectively. One can see that the profiles develop oscillations as (3 increases, 
but that the behaviour close to the bottom (in log scale) remains the same. In panels (d) and 
(e), we plot the coefficients A, B, C and D vs (3 (still for rjo = 10 -4 ). They are weakly affected, 
meaning again that the behaviour close to the bottom is almost unchanged. 



5.1. A second order turbulent closure 

As already stated, to solve quantitatively the 'dune problem', we need to take correctly 
into account the effects inducing a phase shift between the stresses and the relief. As a 
matter of fact, a first order closure assumes that the turbulent energy adapts instanta- 
neously to the mean strain tensor. To take into account the lag between the stress and 
the strain tensors, one needs to formulate a second order turbulent closure. This can be 
achieved by deriving dynamical equations for r^, which relax the stresses towards their 



steady state expression prescribed by equation (2.6) (see Appendix H3 



DtUk = d t T ik + UjdjTik = 



l7l 



k 2 L 2 



nnij 



(5.1) 



In this relation, the parameter (3 encodes the time lag between an increase of the mean 
shear strain rate and the corresponding re-adjusment of the fluctuations of the shear 
stress. We expect (3 to be on the order of unity. 

In figure [9] we show the effect of this new parameter. As expected for inertial effects in 
a relaxation process, finite values of (3 generate oscillations in the vertical profiles of the 
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Figure 10. Effect of the normal stress anisotropy. (a-d) Coefficients A, B, C and D as a 
function of r/o for different values of \x ~ x\ (0.1, 1, 5 and 10). The dashed lines correspond to 
the isotropic case. Inset (e): ratio B/A. Arrows indicate increasing normal stress anisotropy. 



velocities and stresses. The example of the horizontal velocity is displayed in the panels 
(a), (b) and (c). The amplitude and the frequency of these oscillations increase with 
(3. Interestingly, these oscillations do not affect much the behaviour of the modes close 
to the bottom. As a consequence, the coefficients A, B, C and D are weakly affected 
by (3, see panels (d) and (e). Interestingly, both A and B decrease as /3 increases and 
their ratio remains roughly constant. (3 has thus a negligible effect on the emergence of 
bedform, and we shall keep it to zero in the rest of the paper, as well as in part 2. 



5.2. Reynolds stress anisotropy 

It is an experimental fact that, in a turbulent boundary layer close to a rough wall, the 
Reynolds stress tensor is not isotropic: t xx is significantly larger than the other compo- 
nents ( |Raupach et al. 1991[IShafi fc Antonia 1995| ). Besides, anisotropy seems less pro- 
nounced for a larger bottom roughness ( |Krogstad fc Antonia 19 94 , Kcirs bulck et al. 2 002 ) , 
an issue which is however still matter of debate (Krogsta d et al. 2005| ). 

To account for this Reynolds stress anisotropy, it is easy to generalise the Prandtl-like 
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first order turbulent closure (2.6 1 with the following expression: 

'1 



= K 2 L 2 \j\[^\j\S.. 



Hi 



(5.2) 



where the value of \x now differs from that of \ z . Following the above-cited literature, 
we expect Xx/xi to be around 1.3-1.5. The modification of the matrix V due to this new 
closure is detailed in Appendix [A] It is shown that the relevant anisotropic parameter 
entering the equations is x\ ~ xl i f° r which a realistic value is on the order of unity. As 
evidenced in figure [10] the corresponding values of the functions A, B, C and D are not 
much affected by this anisotropy in the relevant range of r/Q. This is particularly true for 
the coefficients C and D, as well as for the ratio B/A as soon as r/o < 10 -2 . The normal 
stress anisotropy has thus a negligible influence on ripple and dune formation, and we 
will not take it into account in the rest of the paper, as well as in part 2. 

5.3. A moving bottom 

In order to investigate the effect of a moving bottom on the stress coefficients, we consider 



a bottom profile of wavevector k like in (3.1 ), but which is now function of both position 
x and time t: 

Z(x,t) = (e at e l{kx ~ UJt K (5.3) 

In this expression, a represents the growth rate of the profile, and uj/k its phase velocity. 
As discussed in Colombi ni fc Stocchino 2005} this investigation is important as we wish to 
use the present hydrodynamical study in the context of the formation and development 
of bedforms, which do have a (small) growth rate and a (small) velocity. Following 
expression (5.3 1, we modify those for t he funct ions U, W, St and S n by inserting the 
extra-factor e^ -1 ")*, along the lines of ( 3.3|3.6 1. 

In this new case, the linearised equations (3.171 and (3.181 of section [3] must then be 
modified in the following manner: 



S' t = 



a 



- i- h ifi ) U + f/W + iS n + iS x 



iS z 



S, n — 



o~ 

kit* 



-i- )W + iS t . 



(5.4) 
(5.5) 



The linear equation (3.241 is then the same, but now with the modified matrix 







ku* 



%\n- 
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(5.6) 



With the surface layer model described in section |4.1| the non-slip boundary conditions 
on the bottom can be written as 

a . lo 



U(0) = -m'(0) 



ct 



W(0) = 



Jen* 



(5.7) 



The result of the integration of this new system is shown in figure [TT] for the shear 
coefficients A and B. One can see that departure from the static case a — and to = is 
noticeable only for values of if— and of order one. The effect of the wave propagation 
of the bedform can be understood by a simple argument. When the bedforms propagate 
upstream (to < 0) the relative flow velocity seen by the structure is larger so that it 
induces a larger shear stress modulation. As A + iB is by definition the basal shear stress 
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Figure 11. Effect of the motion of the bottom. Stress coefficients A and B as a function of r/o 
for different values of the bottom growth rate a (panels (a) and (b)), and different values of the 
bottom pulsation u) (panels (c) and (d)). Arrows indicate increasing values of a and u. In panels 
(a) and (b), grey dotted lines correspond to = 10, 5, 2, 1 and 0.5, the black dotted line 
being for ^S- — 0.1. The black dashed line corresponds to r^— = —0.1, the grey dashed lines 
being for rf— — —0.5, —1, —2, —5 and —10. In panels (c) and (d) grey dotted lines correspond 



to 
to 



2 and 1, the black dotted line being for 



0.1. The black dashed line correspond 



-0.1, the grey dashed line being for -rr— = —1 and —2. For comparison, in all panels 



the solid lines correspond to the static case a — 0, uj = 0. 



rescaled by v%, both A and B get larger. Reciprocally, when the bedforms propagate 
downstream (w > 0), these coefficients are reduced. Consistently with this argument, 
the ratio B/A is only weakly affected by Li (not shown). The growth rate a affects A 
and B in opposite ways and thus changes the phase shift between the shear stress and 
the topography. For a > 0, A is increased while B is reduced. We have not been able to 
interpret this behaviour in a simple way. 

As discussed in part 2, for ripples in water flows these dimensionless numbers are 
respectively on the order of 10 -3 and 10 -2 . They would be even smaller for bedforms of 
larger wavelength. As a consequence, the assumption that the bottom can be treated as 



fixed for the study of bedforms is definitively valid (see also the discussion of figure 22 
below) . 
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6. Weakly non- linear expansion 

In this section, we investigate the non-linear effects at finite values of the rescaled 
bottom corrugation amplitude In particular, we wish to address two issues: we wish 
to relate the hydrodynamic roughness to geometrical quantities; we aim also to describe 
the first non- linear corrections to the basal stress coefficients A, B, C and D. These 
results will be used in the second part of this work, to determine the equilibrium height 
of ripples. In this context, few authors have studied the influence the non-linear terms 
from hydrodynamics on the shape flFreds0e 1974ft or the stability ([Ji & Mcndoza 1997) 
of the bedforms. In contrast to the present paper, however, both of these works de- 
scribe the turbulent closure with a constant eddy viscosity. In a very empirical manner, 
Mc Lean fc Smith 19861 computed the flow over two-dimensional bedforms of arbitrary 
amplitude by the use of a wake function, as described in Schlichting fc Gersten 2000[ 
coupled with a potential flow modified to take into account the velocity logarithmic law. 
Finally, the linear results of Colombini 2004 have also recently been extended to the 
weakly non-linear situation (Colombini & Stocchino 2008). 

6.1. Expansion in amplitude 

For this purpose, we perform an expansion with respect to the bottom corrugation am- 
plitude, and introduce non-dimensional the following functions for the different orders: 

u x =u*[ti+ {kCy^Ut + (k() 2 U + (k0 2 e 2ikx U 2 + (k() 3 e* kx U 3 ] , (6.1) 
u z = u, [(k(y kx W! + (k() 2 W Q + (kQ 2 e 2lkx W 2 + (kQ 3 e ikx W 3 ] , (6.2) 
r xz = -ul [1 + {K)e lkx S tl + (kC) 2 S t0 + (k() 2 e 2lkx S t2 + (k() 3 e tkx S t3 ] , (6.3) 
P + r zz = Po + ul [(k()e lkx S nl + (k() 2 S n0 + (k() 2 e 2lkx S n2 + (k() 3 e lkx S n3 ] , (6.4) 
r zz - r xx = ul [(k(y kx S dl + (k() 2 S d0 + (K) 2 e 2lkx S d2 + (k() 3 e lkx S d3 ] . (6.5) 

Terms in (k() 3 e 3lkx , which do not contribute to the harmonic response (i.e. to the stress 
coefficients), as well as terms of higher order than (fc£) 3 are neglected. Although the 
principle of the expansion in amplitude is simple, the actual calculations are painful, 
and the technical details have been gathered in appendix [E] In summary, the non-linear 
effects result from the expansion of the mixing length (terms in (&C) 2 ) an d from the self- 
interaction of the linear perturbations: in particular, the combination of terms (k(,)e lkx 
generates second order terms in {kQ) 2 . All involved functions are complex, except /x and 
those related to the second order homogeneous corrections (index 0). To avoid the deter- 
mination of the asymptotic behaviours in this case, we have chosen to treat the surface 
layer by the simple phenomenological model of section [4] Plugging the above expressions 
into the Navier-Stokes and turbulent closure equations, one eventually obtains a linear 
hierarchy of linear differential equations: 

^X a = V a X a + S a , (6.6) 
ar] 

where X a = (U a ,W a , St a , S na ). Of course, V\ and S\ are the matrix and vector of 



expression (3.241. We find that V 3 =V\ 1 whereas V 2 is slightly different and Va is very 
simple: 

/ -i 0\ 



Vi 



-i 
m+^r) pt i 



-ni % J 



(6.7) 
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Figure 12. Vertical profiles of the homogeneous second order corrections Uq (a) and Sto (b). 
These curves have been computed with r/o = 2. 10 -3 . Uo tends towards a negative constant 
value —E at large £, which corresponds to an increased roughness at large distance from the 
wall. Notice also that, close to the bottom, Sto has a constant value, reminiscent of the inner 
layer. 
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(6.8) 



All the heaviness of the method is encoded in the expressions of the right hand terms S a : 
the components of such vectors at a given order depend on the lower order functions X a 
and their derivatives. The integration has thus to follow the hierarchy of the equations, 
one order after the other. 

6.2. Boundary conditions 

The boundary conditions must be specified in order to perform the integration. Following 
the geometrically induced roughness model (section we require that both components 
of the velocity should vanish on the bottom. These conditions express easily in the shifted 
representation, i.e. written in terms of the curvilinear coordinates (see appendix[C|: they 
simply read U a (0) = and W^O) = 0. In terms of the Cartesian functions, we get: 



Ui(0) 
U o (0) 

U 2 {0) 

u 3 (o) 



A*'(0), 

V (0) - ^ M ' 2 (0) - i//(0) [S tl (0) + s t \(o)] 



1 



I M '"(o) + ! M '(o)-^y 3 (o) 



(6.9) 
(6.10) 

(6.11) 



1 
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[/i"(0) + 2 KA / 2 (0)] [2S tl (0) + S&(0)] 
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/x'(0)5 tl (Q) [5*1(0) + 2S; 1 (0)} - ^M'(O) [2S nl (0) + S* nl (0)} 
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Figure 13. (a) Vertical profiles of the homogeneous component of the velocity fi + (k£) 2 Uo for 
fc£ = 0, kC = 0.2 and k( = 0.3. The extrapolation to of the upper part of the curves gives the 
value of the hydrodynamical roughness seen at a distance from the bottom, (b) Coefficient E 
as a function of 770. A good fit of this function is given by E — 2.75 (In 770 — 0.62) 2 . 



V(0) [2S«,(0) + S t2 (0)} , 



(6.12) 



and 



W x (0) = 0, 
W (Q) = 0, 

W 2 (0) = -- 
W 3 (0) 



m'(o), 



i M "(o) - |^' 2 (o) - ^M'(O) [25 tl (0) + SM0)] 



(6.13) 
(6.14) 

(6.15) 
(6.16) 



As in previous sections, for each order, the solution is expressed as a linear superposi- 
tion of the form: X a — X sa + at a X ta + a na X nai where the different vectors are solutions 
of the following equations: 



drj 



V a X s 



S a 



drj 



-X t , 



drj 



-X„ 



v a x t , 



*PaX n 



with 



with 



with 



X sa (0) 



x ta (o) 



X na (0) 



( u a (o) \ 

W a (0) 


V j 

(°\ 



1 

f°\ 




V 1 / 



(6.17) 



(6.18) 



(6.19) 



Again, for the top boundary conditions, we introduce a lid at finite height H, impose 
Staivn) = and W a (r] H ) — 0, and look at the limit H — > +00, i.e. when the results 




Figure 14. Vertical profiles of the third order corrections to the stresses: St3 (a) and S n 3 (b). 
The solid lines represent the real parts of the functions, whereas the dashed lines represent the 
imaginary ones. These curves have been computed with rjo = 2. 10 -3 . Notice again the inner 
layer where the stresses are constant. We write St3(0) = A3 + 1B3 and 5*n3(0) = C3 +1D3. Note 
that both A3 and B3 are negative. 



become independent of H. These conditions lead to two equations on a ta and a na , whose 
solutions give S ta (0) and S na (0) respectively. 

6.3. Results 

We first consider the corrections to the homogeneous base solution (index 0). The corre- 
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The 



sponding velocity profile Uq and shear stress profile Sto are displayed in figure 
term in (&C) 2 of the velocity decreases continuously from z ~ zq to z ~ A and tends 
towards a negative constant — E far from the bottom. Correspondingly, the shear stress 
decreases and tends to far from the ground, as requested. The calculation thus predicts 
an increase of the turbulent drag (i.e. of the basal shear stress) with the corrugation am- 
plitude, due to the non-linearities. In terms of the velocity profile, this translates into 
a hydrodynamic roughness z g of geometrical origin. Identifying the expression of the 
velocity profile far from the bottom with ^ In j-, we simply get: 

In z g = lnz a + K(k() 2 E, (6.20) 

As a consequence, z g increases with E and with the aspect ratio kC,. For the seek of 
illustration, several vertical profiles of the homogeneous part of the velocity /j, + (k() 2 Uo 
are plotted in figure 13 'a) for different values of By definition, z g is the extrapolation 



of the asymptotic part of the curves to vanishing velocities. 

Interestingly, the large scale roughness z g cannot be related to a single geometrical 
length, namely to the corrugation amplitude ( (|Schlichting fc Gersten 2000[|van Rijn 1983 
|Raupach et al. 1991[ [Wibcrg & Nelson 1992). In particular, we predict that the macro- 
scopic roughness z g associated to a wavy surface still depends on the microscopic rough- 



ness zq: as shown in figure 13 ^b), for a given aspect ratio, the apparent roughness z g is 



larger for smaller 770. Furthermore, expression (6.201 is consistent with numerical obser- 
vations of |Taylor et al. 1989| and computations of I Jacobs 19891 who respectively report 
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Figure 15. Third order stress coefficients A3, B3, C3 and D3 as a function 770. Comparing 
the signs with those of A = Ax, B = By, G = C\ and D = D%, it can be inferred that the 
non-linearities oppose the linear effects. In particular, as the amplitude increases, the point of 
maximum shear stress drifts downstream, i.e. B\ + B3 (fcC) 2 decreases. 



linear (in the range 10 -8 < rjo < 10 -3 ) and quadratic variations for the coefficient E as 
a function of In 770. Here, the best fit of this function gives E ~ 2.75 (In 7/0 — 0.62) 2 
The first non-linear corrections to the harmonic terms scale on (fcC) 3 - In figure 



14 



we 

show the corresponding profiles for the stresses. As in the first order case, there exists 
a layer close to the bottom where the stresses are almost constant (inner layer). As 
requested, both components vanish far from the ground. We note St3 (0) = A3 + iB% 
and S n3 (0) = C3 + iD 3 the shear and normal stresses acting on the boundary. These 
coefficients are plotted as a function of 770 in figure [15] Both A 3 and B 3 are negative 
while A = A\ and B = B\ are positive, which means that the first non-linearities oppose 
the linear effects. We will show in part 2 that this is responsible for the selection of the 
height of current ripples. 

The calculation of the non- linear corrections allows to determine the range of amplitude 
£ for which the linear approximation is valid. The representation of the linear solution 
with the fixed system of coordinates (x, z) is valid only when £ is much smaller than zq. 
However, the representation of the same linear solution in the curvilinear coordinates 
(x, z — Z(x)) is valid up to ( of the order of the inner layer thickness £. It is important 
to recall that all the descriptions of the flow equivalent at the linear order (including the 
real solution of the fully non-linear problem) can have very different domains of validity. 

This weakly non-linear computation is illustrated in figure 16 which shows the stream- 
lines for different aspect ratios. It can be seen that the separation of streamlines and the 
subsequent formation a recirculation bubble occurs above an aspect ratio of ~ 1/13, in 
agreement with observations. By comparison, the linear calculation, shown in figure [3j 
leads to the emergence of a recirculation bubble for an aspect ratio of 1/6 i.e. twice 
larger. The non-linear corrections are thus essentials to capture quantitatively the flow 
features. Panel 16 computed for an aspect ratio of 1/8 shows a well-developed recircu- 
lation bubble. The distortion of the separation streamline is not realistic, indicating the 
upper limit of validity of the model. Fortunately, the aspect ratio of ripples and dunes 
is typically smaller than 1/10, which falls into the domain of validity of the calculation. 

Several experiments in flumes or wind tunnels have been performed to measure velocity 
and Reynolds stress profiles over two-dimensional fixed symmetric or asymmetric bed- 
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forms flWiberg fc Nelson 1992"llNelson et al. 1993[lMcLean et al. 1994|lBen"net fc Best 19951 

Coleman et al. 2006, Vcnditti 2007). Closer to our calculations, several experiments 

with a sinusoidal bottom are also reported in the literature (e.g. IZilker et al. 1977[ 

Zilker fc Hanratty 1977[|Abrams fe Hanratty 11)85] IBTickies et al. 1984l|Nakagawa fc Hanratty 2001 



Pogg i et al. 2007D . Direct or large eddy numerical simulations of Navier-Stokes equations 



have also been performed in this geometry (e.g. |de Angelis et al. 1997 Henn fc Sykes 1999[ 



Salv etti et al. 2001]) . For comparison with our model, data from Gong et al. 1996 have 



been chosen. They have been performed in a wind tunnel over sixteen waves with a 
wavelength A — 609.6 mm and a trough-to-crest amplitude 2£ — 96.5 mm, covered with 
a carpet to make them aerodynamically rough. This corresponds to an aspect ratio of 
1/6 much too large to be in the domain where the model is quantitative. Unfortunately, 
we have not found any better data-set for the seek of comparison. The vertical profiles 



of the velocity measured at different locations are shown in figure 17 'a). More precisely, 
these authors have measured the average of the instantaneous velocity modulus i.e. a 
quantity that is always positive and that mixes the average velocity and the fluctuations. 
Although a recirculation bubble is present, these data cannot show it. We compare these 
profiles to those computed at the upper limit of validity of the non-linear calculation 
(k( = 0.3). Yet, the agreement with our computation is fair; in particular, the way 
the flow is accelerated over crests and decelerated in troughs in qualitatively well repro- 
duced. The profiles at A/4 and —A/4 from the crest are close to each other, indicating a 
re-symmetrisation of the flow by non-linearities. Note that the slight difference between 
these profiles is qualitatively reproduced by the model. Looking at the upper part of the 
experimental profiles, one sees that they would extrapolate to around 10 mm while 
the ground roughness is slightly smaller than 1 mm. These two roughness' are particu- 
larly visible on the profile measured on the crest. The model is particularly successful in 
reproducing this feature. 

The non-linear effects on the flow over obstacles are often described in terms of 
boundary layer separation. It has been proposed by [Jensen & Zcman 1985 (see also 



Finniga n et al. 1990 ) that one could still use the linear flow calculation in that case, 



introducing a Active surface enveloping the obstacle and the recirculation bubble down- 
stream of it. As such an envelope creates a Active bump maximum displaced downstream, 
it artiAcially moves the point of maximum shear stress on the bump in the same direc- 
tion flKroy et al. 2002") lAndreotti et al. 2002[1 . Although this trick is of practical use to 
simulate dunes, this envelope technique is not based on any Arm theoretical ground. The 
weakly non-linear calculation performed here is thus of extreme interest to incorporate 
non-linear turbulent effects in dune numerical models in a more controlled way. More 
generally, it can be used in any problem in which a good approximation of the mean flow 
is needed at low calculation cost, including separation. For instance, it may find direct 
applications in the control of turbulence around vehicles. An important limit of such 
Reynolds averaged calculation is that they do not take vortex shedding into account. 



7. Effect of a free surface 

In this section, we investigate the effect of the additional presence of a free surface at a 
finite distance H to the bottom. This situation is relevant to the flow above river dunes 
(see part 2). We follow the outline of the section [3] but staying for easiness with linear 
calculations in two-dimensional situations. 
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Figure 16. Streamlines computed with all non-linear corrections considered here on bumps 
such that fc£ = 0.2, 0.3 and 0.4 from top to bottom. The fluid flows from left to right. Note the 
progressive development of a recirculation bubble for larger acpect ratios. 



7.1. River equilibrium 
In the case of a river inclined at an angle 9 on the horizontal, the shear stress must 
balance gravity. It thus varies linearly as t xz = g(z — H) sin 9 and vanishes at the free 
surface. By definition of the shear velocity u*, we also write t xz = v%{z/H — 1). In the 
context of a mixing length approach to describe turbulence, this length should vanish at 
the free surface. For the sake of simplicity, following the discussion of section [4] we take 
L = (z + Zo) yl — z/H. This choice results in a base flow that is logarithmic as in the 
unbounded situation: 



u x = — In 1 + — , (7.1) 

K V z 0/ 

which is consistent with field and experimental observations. The stress balance equation 
along the z-axis allows to get the pressure, which reads: 

2 

11 / Z \ 

P + t zz = Po + g(H - z) cos 9 = Po + — ±- 1 - — ) . (7.2) 

tan o \ n / 

We define the Froude number as the ratio of the surface velocity u sur f acG to the velocity 
of gravity surface waves in the shallow water approximation: 

fE^E -L= U ^ln(l+ H -)= ilnfl + - > ) V^9. (7.3) 



/gH \/gH k \ zo J k \ zq , 

In the literature, the Froude number is sometime defined as the ratio of the mean velocity 
to the velocity of gravity waves. We will justify this choice in the next paragraph. The 
Froude number of natural sandy rivers lies in general between 0.1 and 0.3 as they flow 
on very small slopes. Larger Froude numbers are reached in flume experiments. 

7.2. Disturbances 

In the same manner as in section [3] we consider now a wavy bottom Z = Qe lkx . We note 
again rj = kz and r/n — kH. We write the first order corrections to the base flow as 

u x = u* [n + kQe ikx U] , (7.4) 
u z =u*k(e lkx W : (7.5) 
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Figure 17. (a) Rescaled profiles of the mean velocity modulus, measured by |Gong et al. 1996 
in the case of a rough wavy bottom of aspect ratio ~ 1/6 (k( = 0.5). U TC f is the free stream 
reference velocity. The symbols correspond different longitudinal locations, as shown in the 
schematics below the data (crests □, troughs ♦ and half ways up- V and down-stream •). (b) 
Velocity profile predicted by the present model with a sinusoidal bottom of aspect ratio 1/10 
(kC = 0.31), plotted with the same symbol code. Dotted line: base velocity profile. 
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where, in accordance with equation (7.1 1, the function /i is defined by the relation (4.2|. 



The free surface is also disturbed by the presence of the non-uniform bottom, and we 
denote H+A(x) the flow depth at the position x. The modified expression for the mixing 
length then reads 



(z + z - Z) 



H + A 

WTa~- 



(7, 



Linearising the free surface profile as A(x) = 5(e , one can expand L to the first order 
as 



kL=( V + T] )Jl 




kC,e 



ikx 
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l 

2m 
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>IH 



(7.9) 



The shear stress closure as well as the Reynolds averaged Navier-Stokes equations can 
be linearised in the same way as before, and we finally get at the first order in kC, a 
system of differential equations which can be written under the following form: 



— X = VX + S + 5S 5 , 



(7.10) 
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7.3. Resolution of the linearised equations 



Again, making use of the linearity of the equations, we seek the solution under the form 



X — Xq + atX t + a n X n 



dq 



SXg, where the vector Xg is solution of equation: 

(°\ 





w 



Xs = VXg + Ss 



with 



X S (0) 



(7.13) 



while X , X t and X n are still solutions of equations (4.3 1-( 4.5 1 . The bottom boundary 
conditions U(0) = — 1/(kj7q) and W(0) — are then automatically satisfied. At the free 
surface, we impose the material nature of the surface, W(rju) — ifJ-(vii)S, and vanishing 
stresses: S t (rjH) — and S n (T]ff) = ^/(vh tan 9). These last three conditions select 

the coefficients at and a n as well as the value of S. Finally note that the analytical 
approximation of the solution close to the bottom in the limit 770 — > is the same as in 
the unbounded case - it does not depend on the position of the upper boundary - and 



expressions ( 3.26 )-( 3.29 1 are thus still correct in the limit H ^> zq. 



7.4. Results 

In order to evidence the role of the free surface, we have plotted the stress coefficients 



A and B as functions of t/q in figure 18 for different values of H/zq. For a large enough 



wave-number k (a small enough wavelength A) , one recovers the plots of the panels (a) 
and (b) of figure [6] independently of H/zq. This means that for a bottom wavelength 
much smaller than the flow depth H (i.e. for subaqueous ripples), the free surface has 
a marginal effect and the results of section [3] apply. For smaller rj , however, the curves 
exhibit a peak, whose position depends on the value of H/zq, followed by a diverging 
behaviour when 770 — * 0. As discussed below, this peak can be ascribed to a resonance 
of standing waves at the free surface, excited by the bottom topography, meaning that 
the proper scale is now H and not Zq. As the ratio X/H is the key parameter separating 
ripples from dunes, we shall turn extensively in part 2 to this point. 

The analysis of velocity profiles for different values of kH gives the following physical 
picture. For kH > 1, as for the unbounded case the flow can be thought of as being 
divided into two regions: an inner layer close to the bottom where it can be described 
by the equilibrium approximation and an outer layer behaving like an inviscid potential 
flow, where the profiles can be decomposed into the sum of decreasing and increasing 
exponentials e ±v . For smaller values of kH, this outer region progressively vanishes and 
the whole flow is controlled by the inner layer. 

We display the phase and amplitude of the free surface as a function of kH in fig- 





Figure 18. A and B as functions of r/o for T = 0.9 and H/zo = 10 3 (solid line) or H/zo = 10 4 
(dotted line). In the right part of the plots, the curves collapse on the shape displayed in panels 
(a) and (b) of figure [6] They differ at smaller rjo, showing a resonance peak and a divergence at 
770 -> 0. 



ure 20 T>c). The peak in amplitude accompanied by the phase shift of it are the signature 
of a surface wave resonance. The source of disturbances is of course the corrugation of 
the bed. For kH larger than its resonant value, the bottom and the free surface are in 
phase; conversely, for kH below the resonance, they are in antiphase. In between, at the 
resonance, the phase shift is ip = n/2 (figure [20|( a)) so that the streamlines are squeezed 
downstream to the crest. This resonance is model-independent as it comes from a very 
robust physical mechanism. As the fluid flows over the periodic bottom, gravity surface 
waves are excited at the wavelength A. The latter propagate at the velocity: 

c ~ w surfacc ± J| tanh(fcff) = Ty/^H ± yj | tanh(fcff) (7.14) 

with respect to the bottom (see the friction force model derived in Appendix [F|. As in 
the sound barrier phenomenon, the wave energy induced by the bottom disturbances 
accumulate when this velocity vanishes i.e. for: 



,-^^1 (7.15) 

In the shallow water approximation (kH < 1), this resonant condition gives T = 1 as 
standardly obtained in hydraulics. In the deep water approximation, it gives T = 1/ VkH 
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Figure 19. Relative importance of the inner (white) and outer (grey) layers for kH — 1.65 
(a)-(b), and kH — 0.03 (c), with H/zq = 10 4 . The velocity profiles (bold lines) are compared to 
their asymptotic behaviour in the inner and outer layers (thin lines). The solid lines represent 
the real part of the functions, and the dotted lines the imaginary ones, (a) and (b) show the very 
same profile, but with a logarithmic scale in (a) to emphasize the inner region. The thin lines 
in (a) and (c) represent the asymptotic behaviour in the inner layer. Those in (b) correspond 
to a sum of an increasing and a decreasing exponential of the form exp(±ry), as for an inviscid 
potential flow. 




Figure 20. (a) Streamlines of a flow over a sinusoidal bottom close to the free surface resonance 
conditions {tp — tt/2). The flow is from left to right. Note the squeezing of the lines downstream 
the crest of the bump. Amplitude \S\ (b) and phase <p(6) (c) of the free surface as a function 
of kH for T — 0.8. The peak in amplitude and the phase shift from to n correspond to the 
resonance. The two schematics illustrate the situations in phase or in antiphase. 
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Figure 21. The phase (a) and amplitude (b) of the rescaled free surface deformation 5 = A/£ 
as a function of kH for T — -> (dotted dashed line), T = 0.2 (dashed line), T = 0.4 (thin solid 
line), T = 0.6 (long dashed line), T = 0.8 (dotted line) and T = 1 (solid line), and H/zo = 10 3 . 
Crossing the resonance, the phase shifts from to n. (c) and (d), same for the friction force 
analytical model. 



or equivalently kH = 1/T 2 ( |Kennedy 1963| ). So, for a bottom of wavelength A, the 
flow is subcritical at low T and low kH and supercritical at large T and kH. Ignoring 
dissipation, the Bernoulli relation states that the sum of the gravitational potential energy 
pgA and the kinetic energy 5P u surface ^ s constant along the free surface. The subcritical 
regime corresponds to deep slow flows dominated by gravity: as the velocity increases 
over a bump, the corresponding increase of kinetic energy must be balanced by a loss 
of gravitational potential energy. As a consequence, the free surface is pinched over 
the bump (tp 



7t, see figure 20 



The supercritical regime corresponds to thin rapid 
flows dominated by kinetic energy. By conservation of the flow rate, a pinch of the free 
surface would lead to an increased velocity. As the bump pushes up the free surface, the 
corresponding gain of potential energy should be balanced by a decrease of kinetic energy 
which is achieved by a deformation of the free surface in phase with the bump (ip = 0, see 
figure 20:). In summary, the free surface responds in phase with the excitation at small 
wavelength and becomes delayed as X/H increases. As in a standard second order linear 
system, the disturbance and the system response are in quadrature at the resonance. 
The phase (f and the rescaled amplitude \S\ of the free surface are displayed in fig- 
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Figure 22. Amplitude of the free surface \8\ as a function of kH for different values of the 
growth rate a (a) and the pulsation uj (b) of the bottom boundary. Arrows indicate increasing 
values of a and u>. These graphs have been computed with H/zo — 10 4 and T = 0.8. In panel 
(a), the grey dashed line is for = —1, the back dashed line for —0.1, the black dotted line 
for 0.1 and the two grey dotted lines for 1 and 2. In panel (b), the two grey dashed lines is for 
j^- — —5 and —2, the black dashed line for —0.1, the black dotted line for 0.1 and the three 
grey dotted lines for 1, 2 and 5. For comparison, in all panels the solid lines correspond to the 
fixed case a = 0, u = 0. 



ure 21'a)-(b) for different values of T. In panels 21'c)-(d), they are compared to the 
analytical predictions obtained using a much simpler closure (Appendix [F| . One can see 
that the amplitude of the resonance increases with the Froude number. For very small 
T ', the phase curve is more complicated to interpret, but note that this corresponds to 
a vanishing amplitude 5: the resonance essentially disappears. For kH — > 0, the free 
surface amplitude seems to converge to some finite value, but the phase slowly goes back 
to 0. This gentle crossover is indeed expected at very large wavelength, a situation for 
which the free surface must follow the bottom topography. 

It is interesting to investigate how the resonance is affected by the fact that the bottom 
moves i.e. can grow or propagate. Following the notations introduced in|5.3| we display 



in figure 22 the amplitude of the free surface \S\ as a function of kH for different values 
of the growth rate a and the pulsation u>. For positive growth rates a, the Q-factor 
of the resonance gets smaller but the resonant wavenumber is not affected. A bottom 
propagating at the velocity uj/k moves the resonant peak along the fc-ff-axis. For a 
positive propagation velocity uj/k the surface velocity with respect to the bottom and 
thus the effective Froude number get reduced. As a consequence, the resonant wavelength 
gets smaller - and kH larger. Conversely, the peak moves to smaller wavenumber for an 
upstream moving bottom. As in sub-section |5.3[ these effects are noticeable for values of 



r 2 - and jf— of order one, while realistic values are respectively on the order of 10 2 and 
1CP 3 . For the ripples and dunes problem (part 2), the bedform motion can thus be safely 
ignored in the hydrodynamical calculation and the effect of free surface interpreted in 
terms of surface standing waves. 

The basal shear stress and pressure and subsequently the coefficients A, B, C and D 
are modified by the presence of the free surface when kH is of order one and below. 
In figure |23j the coefficients are plotted as functions of kH for different values of the 
Froude number. One can see that the resonance peak is more pronounced for larger T 
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- they are actually not visible when T is too small. In agreement with the streamlines 
of figure |20[ a), which shows a squeezing downstream the bump crest, the peak of B 
is negative, corresponding to a phase delay of the stress with respect to the bottom. 
Furthermore, the curves corresponding to the presence of a rigid lid at the same height 
H do not exhibit these peaks. Finally, the diverging behaviour of B as kH is also 
a free surface effect as, in the same limit, B reaches a plateau in the case of a rigid top 
boundary. The behaviours of A, B, C and D at small kH (below the resonant condition) 
can be determined analytically using the simple closure proposed in the Appendix [F] For 
small kH, we get A oc l/(kH), B oc -\/{kH) 2 , C oc l/(kH) 3 and D oc l/(kH) 2 . These 
scalings fit fairly well the solutions of the full equations. 

As a conclusion, there are two situations in which the excitation of standing waves by 
the topography affects significantly the characteristics of the inner layer: (i) around the 
resonance, since the surface wave amplitude is very large and (ii) for vanishing kH, when 
the distance H between the topography and the free surface becomes so small that the 
inner layer invades the whole flow. 

8. A qualitative summary of the results 

As this article is based on a rather technical ground, it is useful to sum up, in a 
qualitative manner, our key results and to put them in perspective with respect to the 
second part of the paper. In the context of the formation of ripples and dunes from a 
flat sand bed submitted to a turbulent flow, a central issue is the description of the basal 
shear stress, which controls bed load transport. Due to the scale separation between the 
typical evolution time of the bedforms and that of the flow, the bottom can be considered 
as quasi-static. The hydrodynamics can be investigated independently of the transport 
issue. 

In the traces of the seminal work of Jackson fc Hunt 1975[ three main regions can be 
evidenced in the turbulent flow over a wavy bottom. 

• An outer layer, away from the bottom, in which the flow is well described by inviscid 
potential equations, i.e. where the perfect flow approximation is valid. The streamlines 
follow the topography so that the velocity is in phase with the bottom. 

• An inner layer, which corresponds to the region where the inertial terms of the 
Navier-Stokes equation are negligible, and thus where the longitudinal pressure gradient 
is balanced by the transverse mixing of momentum due to turbulent fluctuations i.e. 
by the Reynolds shear stress transverse gradient. The thickness t of the inner layer is 
related to wavelength by A ~ i ln 2 (^/zo). At the transition between the inner and outer 
layers, the fluid velocity is slowed down by the shear stress. Due to inertia, the velocity 
is always phase delayed with respect to the shear stress. As the velocity is inherited from 
the outer layer, the shear stress is phase-advanced with respect to the topography. 

• A thin surface layer of thickness h , which is responsible for the hydrodynamical 
roughness zq seen from the inner layer. The dominant physical mechanism at work in 
this surface layer can be of different nature. For instance, zq can result from the mixing 
due to roughness elements, the predominance of viscous dissipation, or the presence of 
bed- load transport. Specific scaling laws for zq and ho are obtained for each of theses 
cases. 

The linear relationship between the stresses and the bottom profile can be encoded 
into two complex coefficients A + iB and C + iD. The tangent of the phase between 
the stresses and the bottom is given by the ratios BJA and D/C respectively. These 
coefficients are key inputs for the sediment transport issue (see part 2), and are needed 
to compute the dispersion relation of the bedforms, i.e. the growth rate as a function of 
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Figure 23. A, B, C and D as functions of kH, for T = 0.1 (dotted dashed line), T = 0.8 
(dotted line) and T = 1 (solid line). The dashed lines correspond to a rigid boundary at the 
same height H. The plots have been computed for H/zo = 10 3 . Comparing a free surface to a 
rigid boundary condition (or to the case H <C A), it can be inferred that the hydrodynamics is 
controlled by the surface waves. In particular, the resonance leads to a drop of the shear stress 
component B i.e. to a downstream shift of the point of maximum shear stress.) 



the bedform wavenumber, which tells whether a given wavelength is stable or unstable. 
Their dependencies with respect to several important parameters such as the length scale 
ratios Zq/X or H/X have been determined. Their sensitivity to different ways of modelling 
turbulence or of imposing the bottom boundary conditions have also been discussed. Our 
conclusions regarding these functions can be summed up as follows: 

• A and B are generically positive, whereas C is positive and D negative. This means 
that the shear stress profile reaches its maximum before the crests of the bumps. This 
effect can be visualised on the streamlines, which are squeezed at this pint of maximum 
shear. For the normal stresses, the situation is opposite: the pressure maximum is slightly 
delayed after the crests. As the pressure is almost constant across the inner boundary 
layer, this phase shift D/C is less pronounced by an order of magnitude. A, B, C and D 
have weak dependencies on the ratio X/z . The shear stress phase shift B/A vanishes for 
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asymptotically small kzo and gently increases with ln(fczo). As the inner layer thickness 
£ becomes of the order of the surface layer thickness h , the phase shift drops. 

• These features are robust to (i) the turbulent closure, (ii) the existence of a Reynolds 
stress anisotropy, (iii) the motion (growth or propagation) of the bottom. 

• Much more important is the role of the free surface in the case of a water depth H 
on the order of the bottom wavelength A. The undulations of the bottom excite standing 
gravity waves at the free surface. Resonant conditions are reached when these waves pre- 
cisely propagate at a velocity equal to that of the flow, i.e. when T 2 ~ tanh(fciJ) / (kH) . 
At the resonance, the response of the free surface is in quadrature with the disturbance. 
When the Froude number is large enough, the deformation \6\ of the free surface at the 
resonance is so large that it has a strong effect on the flow close to the bottom. In 
particular, the streamlines are squeezed downstream the crests of the bump, so that the 
shear stress becomes phase-delayed with respect to the topography (B < 0). 

• Another effect due to the presence of a free surface is found for kH — > 0. In this limit, 
the water depth becomes very thin in comparison to the wavelength, and the inner layer 
invades the whole flow (i.e. t ~ H). In this situation, the shear stress and the bottom 
profiles tend to be in phase so that B/A — > 0. Moreover, the shear stress becomes phase 
delayed (B < 0) below a threshold value of kH that increases with the Froude number. 

• The shear stress profile is insensitive to the mechanisms at work in the surface 
layer provided that its thickness ho is smaller than the inner layer thickness t: the 
hydrodynamical roughness Zq is the single quantity inherited from the surface layer. 
The asymptotic calculation performed by IJackson fc Hunt 19751 is recovered but only 
for asymptotically large ln(A/zo)j & limit hardly reached in real problems. When ho 
is comparable or larger than £, the shear stress coefficients A and B are smaller in 
the hydraulically smooth regime where the roughness is due to viscosity than in the 
hydraulically rough regime. In the later case, the mixing of momentum in the surface 
layer is dominated by the turbulent fluctuations induced by the roughness elements. In 
the ripples and dunes problem, these results directly apply to the case where sediments 
are transported with negligible feedback on the flow (erosion limited transport, see part 
2). When the extra-stress due to bedload transport is significant (momentum limited 
transport), we predict that A and B become larger, but with an almost identical phase 
shift B/A. 

• The normal stress profile is almost independent of the surface layer model. 
Beyond the linear case, we have expended the hydrodynamical calculation to the third 

order in In contrast to the linear calculation, the separation of streamlines, when 
the shear stress and the pressure gradient are antagonist, and the associated formation 
of a recirculation bubble are obtained for realistic values of the aspect ratio. At the 
quadratic order, we get a correction to the mean velocity profile which corresponds to 
a roughness z g at large scale due to bottom corrugation. This roughness of geometrical 
origin depends on the 'microscopic' roughness zo inherited from the surface layer and 
is not simply proportional to the amplitude of the bottom profile. We get the scaling 
law: z g = zq exp(n (k() 2 E), where the factor E weakly depends on X/zq. Looking at 
the shear stress modulation, we have shown that the phase shift between the shear stress 
and the bottom profile is reduced by non-linear corrections and change sign at some 
particular aspect ratio. In the second part of this paper, this non-linear analysis is used 
to predict the selection of the aspect ratio of mature ripples and dunes. 



The understanding of the qualitative reason for the upstream shift of the maximum 
shear stress on a bump has been meditated with A.B. Murray. This work has benefited 
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Appendix A. An anisotropic turbulent closure 

The calculations of section [3] can be generalised in the case of the following anisotropic 
stress-strain relationship: 



At the linear order, the velocity, pressure and stress fields read: 

u x = [fi + kCe ikx U] , 
u z = u t kQe ikx W % 



' xz — 1 zx 



-ul [1 + k(e lkx S t ] , 



Po + u- 
1 



X 2 Z + Ke lkx S n 



Xi + k(e l « x S z 



xl + kQe ikx S x 



and the stress equations can be simplified into 



n'S t = 2(U' + iW) - 2k 2 ( v + r) )n' 3 , 
H'S XX = -2iU + 2 -xl{U' + iW) - \xW\ 

»'S ZZ = -2W' + 2 -xl(U' + iW) - \xW 2 - 



The normal stress difference is this time: 

_ -UU 2 X l~xl 

*->XX ^ZZ / "To / 

fi' 3 fi' 
so that on gets the following four closed equations 

1 



(V +iW- Kfi' 2 ) , 



-iW+^'S t + Kfl' 2 , 



U' 

W = -iU, 

S' t =(in+^ju + »'W + l -(xl- xl)St + iS n , 
S' n = -inW + iSt. 



(Al) 

(A 2) 
(A3) 
(A 4) 

(A 5) 
(A 6) 
(A 7) 

(A 8) 
(A 9) 

(A 10) 
(All) 

(A 12) 
(A 13) 
(A 14) 
(A 15) 



As before, they can be written in the usual compact matrix form (3.24), now with 

( ° 



V 















Ifl 



\ 



vt lixl-xl) * 

— /xi i J 



(A 16) 
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Appendix B. A second order turbulent closure 

B.l. Relaxation equation 

The dynamical equations governing the second-order moments Tjfc can be derived rigor- 
ously. Under the assumption of turbulence isotropy at the dissipative scale, it can be 
written under the form: 

o 

D t T ik = d t T lk + UjdjT lk = -r kj d 



i"t " TijdjUk — dj(j)ik — Kik — T^SikS. 



(Bl) 



e is the dissipation rate; 



u i u j u 'k ^ S ^ ne s P a tial flux of turbulent energy induced 



by fluctuations; the pressure term 7T,fe = u' k dip' + u' i dkP > conserves energy and is usually 
responsible for the isotropisation of fluctuations. 

We wish to get a stress tensor that relaxes towards its steady state expression pre- 
scribed by equation (2.6 1. For dimensional reasons, we write the relaxation rate under 
the form |7|//3, where (3 is a phenomcnological constant, and keep the mixing length L 
fixed by the geometrical distance to the wall. The second moment equation then takes 
the form of a first order relaxation equation: 



D t T ik = d t T ik + UjdjTik 



|7| 




k 2 L 2 



(B2) 



Setting — 0, one recovers the stationary solutions (2.6 1. A finite value of introduces a 
lag between a change of the flow velocity field and the point/time at which the Reynolds 
stress readapts to this change. 

B.2. Equations for 2D steady flows 
For 2D steady situations, the stress relaxation equations are the following: 



u x d x T xz + u z d z T xz = \ \-K 2 L 2 \j\% z 





u x d x T zz + u z d z T z 



M 


ItI 




K 2 L 2 \j\j xx +-K 2 X 2 L 2 \M 



K 2 L 2 \j\j zz + - K 2 x 2 L 2 \j\ 



At linear order, they simplify into: 

(// + i(3fi)S t = 2(U' + iW) - 2k 2 ( V + 77 )^' 3 , 
(// + i0»)S xx = -2*17 + \x 2 {U' + iW) - ^xV 2 , 
(// + t0^)S zz = -2W + ^x 2 (U' + iW) - 2 -x 2 ^ 12 . 



Taking the difference of equations ( B 7 1 and ( B 8 1, one can compute 

-4iU 



to obtain four closed equations: 



U^^iW+^tst + K^, 



W = -iU, 
S' t = [in 4 



fj,' + iftii 



u + n'w + iSn, 



(B3) 
(B4) 
(B5) 

(B6) 
(B7) 

(B8) 
(B9) 

(BIO) 
(Bll) 

(B12) 
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As before, they can be written in the usual compact matrix form (3.24), now with 
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Appendix C. Representation of the disturbances 

C.l. Linear order 

Recall that Z(x) = £e lkx is the bottom profile whose wavelength is A = 2ir/k, and 
rj = kz the dimcnsionlcss vertical coordinate. At the linear order, we write all the 
relevant quantities under the form: 

f = f{r,) + kte ikx f 1 {r } ). (CI) 

An alternative is to use the curvilinear coordinates £ = rj — kZ and write the field / as: 

f = m + kte ikx h{Z). (C2) 

We call these expressions respectively the 'non-shifted' and 'shifted' representations of 
/. They lead to the same linearised equations as they are related to each other at the 
linear order by 

/i = /i+/'. (C3) 
In practice, this is especially relevant for U, for which / — fj, is not constant: U = U+fi' is 
the shifted representation of the first order correction to the horizontal velocity. However, 
W = W, St = St and S n — S n . Importantly, note that the range in 77 for which these 
two representations are valid is not the same a priori. 

C.2. Representations for the non-linear expansion 

All fields are expanded up to the third order in fc£, neglecting also non-harmonic terms 
in (k() 3 e kx . Non-shifted representation of the streamwise velocity: 



Ux 



Uto) + (kOe^U^r,) + (kCyUoin) + (fcC) V lto £/ 2 (7j) + {k^Ye^U^). (C4) 



Shifted representation of the same quantity: 

^ = M (0 + {kOe ikx Ui(0 + (K) 2 Uo(0 + (kC) 2 e 2 ^U 2 (0 + (fcC) 3 e fe f/ 3 (0- (C5) 

Expanding the functions f a (v~kC elkx ) with respect to fc£, one can relate a representation 
to the other as: 



Uo = U + \n"-\(Ui + U?), 



u 2 = u 2 + \li 
1 



(C6) 
(C7) 

(C8) 



1 1 -. 



1 -. 
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Conversely: 

Ux = U 1 +[i', (CIO) 
Uo = U + ~n" + ±{U[ + U[*), (Cll) 

U2 = U 3 + ~n»+±Ui, (C12) 

U 3 = U a + + \u'{ + l -U>r + U' a + \u> 2 . (C 13) 

The passage from a representation to the other for the other fields works the same, except 
that there is no zeroth order (function fi) in the expressions. 



Appendix D. Stream function 

To compute the streamlines, we introduce the so-called stream function ^(x, z), defined 
by d^/dx — —u z and d^/dz — u x . This function is such that u- = 0, so that the iso- 



contours = Cst precisely show the streamlines. Using the continuity equation (2.7 1, 
it is easy to show that a solution is \& = Jdzu x . This integral is computed between 
z = Z (the bottom) and z — z. We note £ = rj — kZ the rescaled distance to the bottom. 
Restricting to the linear order, with the relation U\ — iW[ (equation (3.21 1), we end up 
with 



(Dl) 



where the function /i is given by expression (4.2 1 



In the situation with a free surface, one can use the following representation for the 
field /: 

f = m + Ke lkx m), with ^ = VH H Z +A Z _ Z - (D2) 

This curvilinear variable £ vanishes on the bottom z — Z, and £ = 770 at the surface 
z = H + A. The new function /1 is related to those of the non-shifted representation / 
and /1 as: 

i 1 (D3) 



/i(0 = /i(£)+(i + (*-i)^)/'(0- 



For / = u x , we have / = /1 and fx = U\ = iW[. Consequently, the new stream function 
is 



11 I I 



■/ kx 



+ + (*-!) 



(D4) 



One can check that the free surface is indeed a streamline itself, as one of the top boundary 
conditions is Wi(jjh) = 

In the non- linear (and unbounded) case, we get: 




£ ..1,1 

u (0<% + - t i , + -(u 1 + uz) 



4 r 4 
1 / , 1 



+ (fcC)V te (iw 3 + l»" + \u[ + \u[* + u Q + l -u^j | 
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Appendix E. Weakly non-linear calculations 

Definition of the different functions involved in the expansion: 

u x = u* [/i + (fcC)e te C/i + (k() 2 U + (k() 2 e 2 * kx U 2 + {kQ 3 e ikx U 3 ] , (E 1) 
u z = u. [{k()e ikx Wx + {kQ 2 e 2ikx W 2 + (k() 3 e lkx W 3 ] , (E 2) 

t xz = -ul [1 + (k()e ikx S tl + (k() 2 S t0 + {kQ 2 e 2lkx S t2 + {kQ 3 e lkx S t3 ] , (E3) 
P+r zz =p + ul [(k()e ikx S nl + (k() 2 S rM + (kC) 2 e 2lkx S n2 + (Kfe ikx S n3 ] , (E4) 
r zz - r xx = ul [(k()e ikx S dl + (k() 2 S d0 + (k0 2 e 2tkx S d2 + (k() 3 e ikx S d3 ] . (E 5) 
Expansion of the mixing length: 

K 2 (kL) 2 = ±- 2 ^{kQe ikx + yOC) 2 + y {kQ 2 e 2ikx . 

Expansion of the strain tensor components: 

j xx = 2 ((kCy^ill! + (kQ 2 e 2ikx 2iU 2 + (Kfe lkx iU 3 ) , 

^zz — ~'Jxx: 

7xz = n' + (k()e ikx (U[ + iWi) + (kQ 2 U' + {kQ 2 e 2ikx {U' 2 + 2iW 2 ) 
+ (k() 3 e lkx (U! s + iW 3 ), 

which gives for the strain modulus: 

UiUf 



(E6) 



(E7) 
(E8) 

(E9) 



= ti' + (k()e lkx (U[ + iWi) + (k() 2 
+ (k() 3 e lkx 



U' 



We 



2„2ikx 



U 2 + 2iW 2 



U' 3 + iW 3 -^{U[ + iW l ) + ^ 2 {U[*-iW* 1 ) + 



ul 

(E10) 



and then 
K 2 (kL) 2 \j\ 

+ (k(f 



(fcC)e 



//.-./• 



ft 



12 



(U[ + iW 1 )-2n 



-^(U[ + U[* + iW 1 -iW 1 *) + ^ 2 



+ (k() 2 e 2lkx 

+ (k() 3 e lkx 
1 



^ <w[ + [/;• + 2tw, - iw;) - 4 f 2t/j + t/l(2t/ '*,~ Ul) + ui + mw 2 ) 

4 /i y fi J 



(Ell) 



Expressions of the different functions corresponding to the normal stresses t zz — t xx : 

(E12) 



S dl = % 



(U^Ui* - iW*) - U*(U[ + *Wi)) + 2in(U* - U{), 



(E 13) 
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S d2 = -,U 2 + \ui{U[ + iWi) - UkU u 



Sd3 



Ai Pii 

~^U 2 (U[* - iWt) - 8ikU 2 + t K 2 fi'(2U 1 - Uf) + [/?[/? 



ft- 



(E 14) 
(E15) 
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Expressions of the different S, 
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+ 4zU 2 (U{* - iWfi ~ 8kU 2 + k 2 ij!{2U x - U\) + JkUfU? 
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- ^U{{U l 2 + 2iW 2 ) 
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(E 18) 
\ 

(E .9) 



Appendix F. A friction force closure 

Several of the free surface effects can be recovered within a simple friction force model, 
for which analytical expressions of the linear solution of the flow can be derived. In 
particular, the resonance condition as well as the behaviour of the basal stress coefficients 
A, B, C and D for kH — ► can be found and interpreted. 

F.l. Reference state 

We start from the Navier-Stokes equations for a perfect flow, with a crude additional 
turbulent friction term as an approximation of the stress derivatives: 

d x u x + d z u z = 0, (Fl) 

u x d x u x + u z d z u x = -d x p + gsinO - fl^u x , (F2) 

2u 

u x d x u z + u z d z u z = -d z p- gcos6 -tt—fu z , (F3) 



45 



Physically, the force applied to a fluid particle is directly related to the relative velocity 
with respect to the ground. At an angle 9, the following plug flow is an homogeneous 
solution of the above equations: 



( F4 ) 

u z = 0, (F5) 
p = gcos9(H - z). (F6) 

In order to estimate the value of the friction coefficient, one can make use of the fact 
that typical turbulent velocity vertical profiles are logarithmic. However, as the logarithm 

varies slowly when z is much larger than zq, we write u ~ f^dz u x (z) ~ ^ ^ln ^ — 1 J . 

Identifying the shear stress on the bottom as v% = gH sin 6, we finally get with the 



relation (F4| 



a ~[ajz-i) ■ (F7) 

For H/zq in the range 10 3 -10 4 , we get a typical value for on the order of few 10 -3 . We 
now normalize quantities by u and H and get a single non-dimensional (Froude) number: 

T= r^—„ - (F8) 
V gH cos 9 

F.2. Disturbance 

The starting equations can be linearised around the above reference state. Looking a the 
flow over a corrugated bottom Z(x) — Ce , it is easy to show that the solution is of the 
following form 

u x = u + uk(e tkx [-a + e kz + a_e~ kz ] , (F 9) 

u z = uik(e lkx [a+e kz + a-e~ kz \ , (F 10) 

p = g cos 6(H - z) + u 2 (kH - i2Q)^- [a + e kz - a_e" fez ] , (F 11) 

where a+ and a_ must be determined by the boundary conditions. This exponential 
form is characteristic of potential flows. 

F.3. Boundary conditions 

We require that the velocity normal to the bottom vanish. Following the notations of the 
main part of the paper, we define A such that the free surface is at the altitude H + A. It 
is a material line where the pressure vanishes. The three boundary conditions are then: 

u z {z = 0) = iukCe lkx , (F12) 
u z (z = H) =iu5k(e lkx , (F13) 

where, as before, 5 is defined as A(x) = S(e lkx . The constants a+ and a_, as well as 5 
are thus solutions of 

a+ + a- = 1 , (F 15) 

a + e kH + a_e- kH = 5 , (F 16) 
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from which we get: 



A. Fourriere, P. Claudin and B. Andreotti 

a + e kH - a_e- kH = 



(kH - i2Sl)T 2 ' 



a+ 



(kH - i2Q) tanh kH - 



1 

F 2 



1 



(kH-i2fl) - ^ tanh kH 
(kH - i2Q) tanh kH - ^ 
-7^2 tanh kH 



(F17) 

(F18) 
(F19) 



(kH - i2Q) 

F.4. Basal shear stress and pressure 

The shear stress is not part of the variables of this model, but we can consistently define 
it as t = — flu 2 . Looking at the shear stress Tf, and normal stress pb on the bottom, in 
accordance with the notations of the previous sections of the paper, we introduce the 
coefficients A, B, C and D as 

U = -nil 2 [1 + (A + iB)k(e lkx ] , (F 20) 

p b = gH cos 9 + Vlu 2 (C + iD)k(e lkx , (F 21) 



which gives 



A = 2 



4ft 2 



■pi 



] tanh kH - kH [tanh 2 kH + 1] 



(kH - ^ tanh kH) + 4ft 2 



(F 22) 
(F 23) 

(F 24) 
(F 25) 



It is worth noting that the friction force model predicts negative values of B for any kH. 
This means that there is always a phase delay of the shear stress with respect to the 
bottom, which is a clear disagreement with the full solution. In order to fix this flaw, one 
would need to empirically introduce an imaginary part to ft. Finally, this discrepancy 
shows that a precise description of the phase between the basal friction and the relief is 
a subtle and difficult issue that fully justify the use of a rigorous but heavy formalism. 



B = 


2ft 




[tanh 2 kH - 1] 


't 2 


(kH- 


- ^tanhfci/) 2 + 4ft 2 ' 


C = 


1 


(- 


2ftB\ 


2ft 


kH ) ' 


D = 


1 


(- 


2fL4\ 


2ft 


+ kH • 
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